!Clear SLAE matrix and vector
subroutine clear(ggl,ggu,di,pr,n,m)
    dimension di(n),ggl(m),ggu(m),pr(n)
    real*8 ggl,ggu,di,pr

    call clearVectorDouble(di,n)
    call clearVectorDouble(pr,n)
    call clearVectorDouble(ggl,m)
    call clearVectorDouble(ggu,m)
end
!*************************************************************************
subroutine clearVectorDouble(x,n)
    dimension x(n)
    real*8 x
    do i=1,n
        x(i)=0.d0
    enddo
end
!*************************************************************************
subroutine clearVectorInteger(x,n)
    dimension x(n)
    integer x
    do i=1,n
        x(i)=0.d0
    enddo
end
!*************************************************************************
real*8 function psi1(r1,r2,z1,z2,r,z)
    real*8 r1,r2,z1,z2,r,z
    psi1=((r2-r)/(r2-r1))*((z2-z)/(z2-z1))
end
!*************************************************************************
real*8 function psi2(r1,r2,z1,z2,r,z)
    real*8 r1,r2,z1,z2,r,z
    psi2=((r-r1)/(r2-r1))*((z2-z)/(z2-z1))
end
!*************************************************************************
real*8 function psi3(r1,r2,z1,z2,r,z)
    real*8 r1,r2,z1,z2,r,z
    psi3=((r2-r)/(r2-r1))*((z-z1)/(z2-z1))
end
!*************************************************************************
real*8 function psi4(r1,r2,z1,z2,r,z)
    real*8 r1,r2,z1,z2,r,z
    psi4=((r-r1)/(r2-r1))*((z-z1)/(z2-z1))
end
!*************************************************************************
real*8 function functionX(x1,x2,hx,x,i)
    real*8 x1,x2,hx,x
    k=mu_3d(i)
    if(k.eq.1) then
        functionX=(x2-x)/hx
    else
        functionX=(x-x1)/hx
    endif
end
!*************************************************************************
real*8 function functionY(y1,y2,hy,y,i)
    real*8 y1,y2,hy,y
    k=nu_3d(i)
    if(k.eq.1) then
        functionY=(y2-y)/hy
    else
        functionY=(y-y1)/hy
    endif
end
!*************************************************************************
real*8 function functionZ(z1,z2,hz,z,i)
    real*8 z1,z2,hz,z
    k=ku_3d(i)
    if(k.eq.1) then
        functionZ=(z2-z)/hz
    else
        functionZ=(z-z1)/hz
    endif
end
!*************************************************************************
subroutine calculateValueForNotTermNodes(v,igT,jgT,ggT,n,m,nTerm,nNotTerm)
    dimension v(n),ggT(m),igT(nTerm+1),jgT(m)
    real*8 v,ggT
    do i=1,nTerm
        v(nNotTerm+i) = 0.d0
        do j=igT(i),igT(i+1)-1
            v(nNotTerm+i)=v(nNotTerm+i)+v(jgT(j))*ggT(j)
        enddo
    enddo
end
!*************************************************************************
subroutine mult_sq(g,b,v,n)
    dimension g(n),b(n,n),v(n)
    real*8 g,b,v

    do i=1,n
        g(i)=0.d0
        do j=1,n
            g(i)=g(i)+b(i,j)*v(j)
        enddo
    enddo
end
!*************************************************************************
!a[n]=b[n]+c[n]
subroutine vectorSumDouble(a,b,c,n)
    dimension a(n),b(n),c(n)
    real*8 a,b,c

    do i=1,n
        a(i)=b(i)+c(i)
    enddo
end
!*************************************************************************
!a[n]=b*c[n]
subroutine multVectorByNumberDouble(a,b,c,n)
    dimension a(n),c(n)
    real*8 a,b,c

    do i=1,n
        a(i)=b*c(i)
    enddo
end
!*************************************************************************
integer function searchElement2D(vertex,nVer,r,z,nodesCount,elemsCount)
    integer elemsCount
    real*8 vertex
    real*8 r,z,r1,r2,z1,z2
    dimension vertex(2,nodesCount),nVer(4,elemsCount)

    if(r.lt.1d-3) r=1d-3

    searchElement2D=1
    r1=vertex(1,nVer(1,searchElement2D))
    r2=vertex(1,nVer(2,searchElement2D))
    z1=vertex(2,nVer(1,searchElement2D))
    z2=vertex(2,nVer(3,searchElement2D))
    do while(searchElement2D+1.le.elemsCount.and.(r1.gt.r.or.r2.lt.r.or.z1.gt.z.or.z2.lt.z))
        searchElement2D=searchElement2D+1
        r1=vertex(1,nVer(1,searchElement2D))
        r2=vertex(1,nVer(2,searchElement2D))
        z1=vertex(2,nVer(1,searchElement2D))
        z2=vertex(2,nVer(3,searchElement2D))
    enddo
    if(r1.gt.r.or.r2.lt.r.or.z1.gt.z.or.z2.lt.z)searchElement2D = 0
end
!*************************************************************************
integer function searchElement3D(vertex,nVer,x,y,z,nodesCount,elemsCount)
    integer elemsCount
    real*8 vertex
    real*8 x,y,z,x1,x2,y1,y2,z1,z2
    dimension vertex(3,nodesCount),nVer(14,elemsCount)

    searchElement3D=1
    x1=vertex(1,nVer(1,searchElement3D))
    x2=vertex(1,nVer(2,searchElement3D))
    y1=vertex(2,nVer(1,searchElement3D))
    y2=vertex(2,nVer(3,searchElement3D))
    z1=vertex(3,nVer(1,searchElement3D))
    z2=vertex(3,nVer(5,searchElement3D))
    do while(searchElement3D+1.le.elemsCount.and.(x1.gt.x.or.x2.lt.x.or.&
                                                  y1.gt.y.or.y2.lt.y.or.&
                                                  z1.gt.z.or.z2.lt.z))
        searchElement3D=searchElement3D+1
        x1=vertex(1,nVer(1,searchElement3D))
        x2=vertex(1,nVer(2,searchElement3D))
        y1=vertex(2,nVer(1,searchElement3D))
        y2=vertex(2,nVer(3,searchElement3D))
        z1=vertex(3,nVer(1,searchElement3D))
        z2=vertex(3,nVer(5,searchElement3D))
    enddo
end
!*************************************************************************
subroutine e2d(r,z,vertex2D,nVer2D,vRZ,vXYZ,dimV,nodesCount2D,elemsCount2D,nEl)
    integer elemsCount2D,dimV
    real*8 vertex2D,vertex3D,vRZ,vXYZ
    real*8 r,z,r1,r2,z1,z2,psi1,psi2,psi3,psi4
    dimension vertex2D(2,nodesCount2D),nVer2D(4,elemsCount2D),vRZ(nodesCount2D,dimV),vXYZ(dimV)

    r1=vertex2D(1,nVer2D(1,nEl))
    r2=vertex2D(1,nVer2D(2,nEl))
    z1=vertex2D(2,nVer2D(1,nEl))
    z2=vertex2D(2,nVer2D(3,nEl))

    do j=1,dimV
        vXYZ(j)=psi1(r1,r2,z1,z2,r,z)*vRZ(nVer2D(1,nEl),j)+&
                psi2(r1,r2,z1,z2,r,z)*vRZ(nVer2D(2,nEl),j)+&
                psi3(r1,r2,z1,z2,r,z)*vRZ(nVer2D(3,nEl),j)+&
                psi4(r1,r2,z1,z2,r,z)*vRZ(nVer2D(4,nEl),j)
    enddo
end
!*************************************************************************
subroutine inputParams_task3d(params,sig3d,nMaterials)
    integer sig3d
    real*8 params
    dimension params(10,nMaterials)

    call clearVectorDouble(params,10*nMaterials)
    do while(.not.EOF(sig3d))
        read(sig3d,*)k,(params(j,k),j=1,6),params(10,k)
    enddo
end
!*************************************************************************
subroutine inputParams_task2d(params,sigma_txt,nMaterials2d)
    integer sigma_txt
    real*8 params
    dimension params(*)

    nMaterials2d = 0
    do while(.not.EOF(sigma_txt))
        read(sigma_txt,*)k,params(k)
        if(nMaterials2d<k) nMaterials2d = k
    enddo
end
!*************************************************************************
subroutine getPortraitForEarth(ig,jg,igEarth,jgEarth)
    common /dim3D/nodesCount3D,elemsCount3D,nodesBoundCount3D,edgesTermCount3D,&
    edgesNotTermCount3D,edgesCount3D,nodesTermCount3D,m3D,mTerm3d,nodesNotTermCount3D
    common /dimEarth/nodesNotTermEarth,mEarth
    integer edgesTermCount3D,edgesNotTermCount3D,edgesCount3D,elemsCount3D,elementsEarth
    dimension ig(nodesNotTermCount3D+1),jg(m3D),igEarth(nodesNotTermEarth+1),jgEarth(*)

    k=1
    igEarth(1)=1
    do i=1,nodesNotTermEarth
        j=ig(i)
        do while(j<=ig(i+1)-1.and.jg(j)<=nodesNotTermEarth)
            jgEarth(k)=jg(j)
            j=j+1
            k=k+1
        enddo
        igEarth(i+1)=igEarth(i)+(j-ig(i))
    enddo
    mEarth=igEarth(nodesNotTermEarth+1)-1
end
!*************************************************************************
subroutine rotate(params,nMaterials)
    dimension params(10,nMaterials),mx(3,3),my(3,3),mz(3,3),f(3,3),m(3,3),mt(3,3)
    real*8 params,mx,my,mz,angX,angY,angZ,f,m,mt

    do i=1,nMaterials
        angX=params(4,i)
        angY=params(5,i)
        angZ=params(6,i)
        params(5,i)=params(2,i)
        params(9,i)=params(3,i)
        params(2,i)=0.d0
        params(3,i)=0.d0
        params(4,i)=0.d0
        params(6,i)=0.d0
        params(7,i)=0.d0
        params(8,i)=0.d0

        mx(1,1)=1.d0
        mx(1,2)=0.d0
        mx(1,3)=0.d0
        mx(2,1)=0.d0
        mx(2,2)=dcos(angX)
        mx(2,3)=-dsin(angX)
        mx(3,1)=0.d0
        mx(3,2)=dsin(angX)
        mx(3,3)=dcos(angX)

        my(1,1)=dcos(angY)
        my(1,2)=0.d0
        my(1,3)=dsin(angY)
        my(2,1)=0.d0
        my(2,2)=1.d0
        my(2,3)=0.d0
        my(3,1)=-dsin(angY)
        my(3,2)=0.d0
        my(3,3)=dcos(angY)

        mz(1,1)=dcos(angZ)
        mz(1,2)=-dsin(angZ)
        mz(1,3)=0.d0
        mz(2,1)=dsin(angZ)
        mz(2,2)=dcos(angZ)
        mz(2,3)=0.d0
        mz(3,1)=0.d0
        mz(3,2)=0.d0
        mz(3,3)=1.d0
        f = MATMUL(mx,my)
        m = MATMUL(f,mz)
        mt=transpose(m)

        call matMult(f,mt,params(1,i),3,3,3)
        call matMult(params(1,i),f,m,3,3,3)
    enddo
end
!*************************************************************************
!c[n,k]=a[n,m]*b[m,k]
subroutine matMult(c,a,b,n,m,k)
    dimension a(n,m),b(m,k),c(n,k)
    real*8 a,b,c
    do i=1,n
        do j=1,m
            c(i,j)=0.d0
            do l=1,k
            c(i,j)=c(i,j)+a(i,l)*b(l,j)
            enddo
        enddo
    enddo
end
!*************************************************************************
!*************************************************************************
!*************************************************************************